/*  
 * Contour extraction program
 * 
 * Usage: ./isovis <arg>
 * For the MRI data give an argument, for the CT image give none.
 * Data should be in ../headsq and ../uvadata
 * 
 * Florian Speelman & Jannis Teunissen
 */

#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkVolume16Reader.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkOutlineFilter.h>
#include <vtkCamera.h>
#include <vtkProperty.h>
#include <vtkPolyDataNormals.h>
#include <vtkLODActor.h>
#include <vtkContourFilter.h>
#include <vtkImageData.h>
#include <vtkImageMapToColors.h>
#include <vtkImageActor.h>
#include <vtkImageAccumulate.h>
#include <vtkImageData.h>
#include <vtkXYPlotActor.h>
#include <vtkLookupTable.h>
#include <vtkStripper.h>
#include <vtkCallbackCommand.h>
#include <vtkImageShrink3D.h>
#include <vtkSliderWidget.h>
#include <vtkSliderRepresentation2D.h>
#include <vtkAxisActor2D.h>
#include <vtkDataSetReader.h>

// Bad style can be convenient for simple projects :x
#include "CFList.cpp"
#include "CFMaptors.cpp"

vtkVolume16Reader * readV16Input(void) {
	// Read in the data in the ../headsq directory
	vtkVolume16Reader *v16 = vtkVolume16Reader::New();
	v16->SetDataDimensions (256,256);
	v16->SetImageRange (1,94);
	v16->SetDataByteOrderToLittleEndian();
	v16->SetFilePrefix("../fullHead/headsq");
	v16->SetDataSpacing (0.8, 0.8, 1.5);
	unsigned short mask = (1 << 15) - 1;
	v16->SetDataMask(mask);
	return v16;
}

// Callback for the interaction
class vtkSliderCallback : public vtkCommand
{
	public:
	static vtkSliderCallback *New() 
	{
	return new vtkSliderCallback;
	}
	virtual void Execute(vtkObject *caller, unsigned long, void*)
	{
	vtkSliderWidget *sliderWidget = 
	  	reinterpret_cast<vtkSliderWidget*>(caller);
	if(this->cflist->GetnIso())
		this->cflist->EditIsoValue(static_cast<vtkSliderRepresentation *>(sliderWidget->GetRepresentation())->GetValue(), 0);
	renwin->Render();
	}
	CFList *cflist;
	vtkRenderWindow *renwin;
	
};

int main (int argc, char **argv)
{
	int usingMRI = 0;
	// Read the input
	vtkVolume16Reader *v16;
	vtkDataSetReader *dsReader;
	
	vtkImageShrink3D *shrink = vtkImageShrink3D::New();
	if(argc > 1) {
		usingMRI = 1;
		dsReader = vtkDataSetReader::New();
		dsReader->SetFileName("../uvadata/MRI-Woman.vtk");
		shrink->SetInputConnection(dsReader->GetOutputPort());
	} else {
		v16 = readV16Input();
		shrink->SetInputConnection(v16->GetOutputPort());
	}

	// Shrink the input
	shrink->SetShrinkFactors(3,3,3);

	// Create the renderer, the render window, and the interactor.
	vtkRenderer *aRenderer = vtkRenderer::New();
	vtkRenderWindow *renWin = vtkRenderWindow::New();
	vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
	
	aRenderer->SetViewport(0.0, 0.0, 1.0, 1.0);
	renWin->AddRenderer(aRenderer);
	iren->SetRenderWindow(renWin);
	
	// Set a background color for the renderer and set the size of the
	// render window (expressed in pixels).
	aRenderer->SetBackground(0,0,0);
	renWin->SetSize(1024, 768);

	// Create a class CFList object, to store contour filters
	CFList isoSurfs = CFList(shrink->GetOutputPort(), 10);
	// Set initial isosurface to 500, color red, opacity 0.5
	isoSurfs.AddIsoValue(500, 1, 0, 0, 0.5);
	
	// Create a class CFMaptors object, to store mappers and actors
	CFMaptors maptor = CFMaptors(&isoSurfs, 10);
	for(int i = 0; i < isoSurfs.GetnIso(); i++) {
		aRenderer->AddActor(maptor.GetActor(i));
	}
		   
	// Make a histogram of the data, lots of stuff down here is hardcoded...
	vtkImageAccumulate *accum = vtkImageAccumulate::New();
	if(usingMRI) {
		accum->SetInputConnection(dsReader->GetOutputPort());
	} else {
		accum->SetInputConnection(v16->GetOutputPort());
	}
	accum->SetComponentOrigin(1, 0, 0);
	accum->SetComponentSpacing(30, 0, 0);
	accum->SetComponentExtent(0, 2999, 0, 0, 0, 0);
	accum->IgnoreZeroOn();
	accum->UpdateWholeExtent();
	
	vtkXYPlotActor *plot = vtkXYPlotActor::New();
	vtkAxisActor2D *xaxis = plot->GetXAxisActor2D();
	vtkAxisActor2D *yaxis = plot->GetYAxisActor2D();
	yaxis->SetAxisVisibility(0);
	yaxis->SetTickVisibility(0);
	yaxis->SetLabelVisibility(0);
	yaxis->SetTitleVisibility(0);
	xaxis->SetAxisVisibility(0);
	xaxis->SetTickVisibility(0);
	xaxis->SetLabelVisibility(0);
	xaxis->SetTitleVisibility(0);
	
	plot->SetWidth(.4);
	plot->SetHeight(.3);
	plot->ExchangeAxesOff();
	plot->SetXValuesToValue();
	plot->SetXRange(0, 3000) ;
	plot->AddInput(accum->GetOutput());
	plot->SetPlotColor(0, 1.0, 0.0, 0.0);
	plot->SetPosition(0.6, 0.7);

	aRenderer->AddActor2D(plot);
	
	// Create a slider widget, to adjust the first isovalue
	vtkSliderRepresentation2D *sliderRep = vtkSliderRepresentation2D::New();
 
	sliderRep->SetMinimumValue(0.0);
	sliderRep->SetMaximumValue(3000.0);
	sliderRep->SetValue(500.0);
	sliderRep->GetPoint1Coordinate()->SetCoordinateSystemToDisplay();
	sliderRep->GetPoint1Coordinate()->SetValue(0.765*1024 ,0.75*768);
	sliderRep->GetPoint2Coordinate()->SetCoordinateSystemToDisplay();
	sliderRep->GetPoint2Coordinate()->SetValue(0.965*1024, 0.75*768);
	
	vtkSliderWidget *sliderWidget = vtkSliderWidget::New();
	sliderWidget->SetInteractor(iren);
	sliderWidget->SetDefaultRenderer(aRenderer);
	sliderWidget->SetRepresentation(sliderRep);
	sliderWidget->EnabledOn();

	vtkSliderCallback *callback = vtkSliderCallback::New();
	callback->cflist = &isoSurfs;
	callback->renwin = renWin;

	sliderWidget->AddObserver(vtkCommand::InteractionEvent,callback);
	
	// Create plane datastructures, to later on add a plane to the visualization
	// First define color lookup table
	vtkLookupTable *hueLut = vtkLookupTable::New();
	hueLut->SetTableRange (0, 2000);
	hueLut->SetHueRange (0, 1);
	hueLut->SetSaturationRange (1, 1);
	hueLut->SetValueRange (1, 1);
	hueLut->Build();

	// Create plane structures
	vtkImageMapToColors *planemap = vtkImageMapToColors::New();
	planemap->SetInputConnection(shrink->GetOutputPort());
	planemap->SetLookupTable(hueLut);
	vtkImageActor *plane = vtkImageActor::New();
	plane->SetInput(planemap->GetOutput());

	// Initialize the event loop and then start it.
	iren->Initialize();
	iren->SetDesiredUpdateRate(5.);
	iren->SetStillUpdateRate(0.01);

	// Now enter our command line menu
	double isoValue;
	double rgb[3], opac;
	int menu, n, planeIndex, c[3], firstTime = 1, max[3];
	char dim;

	while(menu != 0) {
		cout << endl << endl;
		cout << "\t# VTK Contour extractor menu\n\n";
		cout << "\t# 0. Exit\n";
		cout << "\t# 1. List current isovalues\n";
		cout << "\t# 2. Add an isovalue\n";
		cout << "\t# 3. Remove an isovalue\n";
		cout << "\t# 4. Go to interactive mode\n";
		cout << "\t# 5. Adjust image quality (performance)\n";
		cout << "\t# 6. Add/edit plane in visualization\n";
		cout << "\t# 7. Remove plane from visualization\n\n";
		
		cin >> menu;
		switch(menu) {
			case 1:
				cout << "\nCurrent isovalues:\n";
				for(int i = 0; i < isoSurfs.GetnIso(); i++) {
					double *c = isoSurfs.GetIsoColor(i);
					cout << "\tContour " << i << " - " << isoSurfs.GetIsoValue(i);
					cout << " RGB: " << c[0] << "," << c[1] << "," << c[2] << " -\tOpac. " << isoSurfs.GetIsoOpacity(i) << endl;
					delete[] c;
				}
				break;
			case 2:
				cout << "Please enter an isovalue\n";
				cin >> isoValue;
				cout << "Enter R,G,B values for the coloring, in range [0,1]\n";
				cin >> rgb[0] >> rgb[1] >> rgb[2];
				cout << "Enter the desired opacity\n";
				cin >> opac;
				isoSurfs.AddIsoValue(isoValue, rgb[0], rgb[1], rgb[2], opac);
				n = isoSurfs.GetnIso() - 1;
				maptor.AddCF(&isoSurfs, n);
				aRenderer->AddActor(maptor.GetActor(n));
				renWin->Render();
				break;
			case 3:
				cout << "\nCurrent isovalues:\n";
				for(int i = 0; i < isoSurfs.GetnIso(); i++) {
					double *c = isoSurfs.GetIsoColor(i);
					cout << "\tContour " << i << " - " << isoSurfs.GetIsoValue(i);
					cout << " RGB: " << c[0] << "," << c[1] << "," << c[2] << " -\tOpac. " << isoSurfs.GetIsoOpacity(i) << endl;
					delete[] c;
				}
				if (isoSurfs.GetnIso() > 0) {
					cout << "Enter the number of the isovalue to remove, -1 to cancel\n";
					cin >> n;
					if( n != -1) {
						aRenderer->RemoveActor(maptor.GetActor(n));
						isoSurfs.DelIsoValue(n);
						maptor.DelCF(n);
						renWin->Render();
					}
				}
				break;
			case 4:
				cout << "Now in interactive mode, press 'e' to exit\n";
				iren->Start();
				break;
			case 5:
				cout << "Give new shrink factors (x,y,z), where 1 means full resolution\n";
				cin >> c[0] >> c[1] >> c[2];
				shrink->SetShrinkFactors(c[0], c[1], c[2]);
				
				max[0] = 256/c[0]-1;
				max[1] = 256/c[1]-1;
				max[2] = 94/c[2]-1;
				switch(dim) {
					case 'x':
					if(planeIndex/c[0] > max[0]) n = max[0];
					else n = planeIndex/c[0];
					plane->SetDisplayExtent(n,n, 0,max[1], 0, max[2]);
					break;
					case 'y':
					if(planeIndex/c[1] > max[1]) n = max[1];
					else n = planeIndex/c[1];
					plane->SetDisplayExtent(0,max[0], n,n, 0, max[2]);
					break;
					case 'z':
					if(planeIndex/c[2] > max[2]) n = max[2];
					else n = planeIndex/c[2];
					plane->SetDisplayExtent(0,max[0], 0,max[1], n, n);
					break;
				}
				renWin->Render();
				break;
				
			case 6:
				cout << "Enter the dimension that the plane is in, as 'x', 'y' or 'z'\n";
				cin >> dim;
				cout << "Enter the index of the plane (" << dim << " = value),\n";
				cout << "the data has dimensions [0,255]x[0,255]x[0,93] so it should be in range\n";
				cin >> planeIndex;

				shrink->GetShrinkFactors(c);
				max[0] = 256/c[0]-1;
				max[1] = 256/c[1]-1;
				max[2] = 94/c[2]-1;
				switch(dim) {
					case 'x':
					if(planeIndex/c[0] > max[0]) n = max[0];
					else n = planeIndex/c[0];
					plane->SetDisplayExtent(n,n, 0,max[1], 0, max[2]);
					break;
					case 'y':
					if(planeIndex/c[1] > max[1]) n = max[1];
					else n = planeIndex/c[1];
					plane->SetDisplayExtent(0,max[0], n,n, 0, max[2]);
					break;
					case 'z':
					if(planeIndex/c[2] > max[2]) n = max[2];
					else n = planeIndex/c[2];
					plane->SetDisplayExtent(0,max[0], 0,max[1], n, n);
					break;
				}
				if(firstTime) {
					aRenderer->AddActor(plane);
					firstTime = 0;
				}
				renWin->Render();
				break;
			case 7:
				if(!firstTime) aRenderer->RemoveActor(plane);
				firstTime = 1;
				renWin->Render();
				break;
			case 8:
				break;
		}
	}
	// If we would not be quitting now, it would have been good to delete the used objects.
	cout << "Bye Bye!\n";
	return 0;
}
